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Abstract 

Biology presents many examples of planar distribution and structural networks having dense sets of closed 
loops. An archetype of this form of network organization is the vasculature of dicotyledonous leaves, which 
showcases a hierarchically-nested architecture containing closed loops at many different levels. Although a 
number of methods have been proposed to measure aspects of the structure of such networks, a robust metric to 
quantify their hierarchical organization is still lacking. We present an algorithmic framework, the hierarchical 
loop decomposition, that allows mapping loopy networks to binary trees, preserving in the connectivity of the 
trees the architecture of the original graph. We apply this framework to investigate computer generated graphs, 
such as artificial models and optimal distribution networks, as well as natural graphs extracted from digitized 
images of dicotyledonous leaves and vasculature of rat cerebral neocortex. We calculate various metrics based 
on the Asymmetry, the cumulative size distribution and the Strahler bifurcation ratios of the corresponding 
trees and discuss the relationship of these quantities to the architectural organization of the original graphs. 
This algorithmic framework decouples the geometric information (exact location of edges and nodes) from the 
metric topology (connectivity and edge weight) and it ultimately allows us to perform a quantitative statistical 
comparison between predictions of theoretical models and naturally occurring loopy graphs. 

1 Introduction 

Among the many different classes of complex systems that can primarily be described as "networks" , an important 
subclass concerns physical networks devoted to transportation of various entities, such as fluids or energy. To some 
extent, structural load-bearing networks can also be considered in this category, as their job is the distribution of 
stress-strain. Besides their evident technological importance, these networks are central to the function of living 
beings; because of their concrete physicality they are sometimes far more accessible to experimental analysis than 
other important biological networks, and hence offer an important window into the organization and function of 
naturally evolved large-scale networks. 

Many biological distribution and structural networks contain dense numbers of reentrant loops. The venation 
of angiosperm leaves (Fig. [T]) [T], the structural veins of insect wings, the continuously adapting foraging networks 
of some fungi and slime molds [2], the vasculature of animal organs such as the adrenal glands, the brain [3] and the 
liver are just a few of a large number of examples where physical networks developed loops in living organisms. These 
networks perform functions crucial to the survival of the organisms that use them. The hierarchical organization 
and the intricacies of the architecture of these highly interconnected networks dictate the efficacy in providing 
support or distributing load under varying conditions. Iir some cases the function of closed loops and how many 




Figure f: Variability in natural loopy networks, (a), (b) Leaf vasculature of two dicotyledonous species, (c) Detail 
of leaf collected from the same plant as leaf (a). The venation of (a) and (c) is predominately reticulate, (b) is percurrent. 
In general, leaves from the same plant (or species) share statistically similar architectural properties, as compared to leaves 
from different species. The scale is f cm. 
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there should be is intuitively obvious; the webbing-like veins of a dragonfly wing have cross-bracings that serve to 
maintain rigidity and resistance with a minimum of weight. In other cases it is not self-evident why there are as 
many loops as observed. 

In many cases, such as leaf venation, loopy networks evolved gradually from a tree architecture [4]. Various 
reasons for the evolution of loopiness in biological distribution networks have been proposed [5l-[7] . These networks 
are the result of developmental processes that frequently dictate not the exact position of each network edge but 
the overall organization in a statistical sense. For example, one can frequently determine by mere inspection of the 
leaf venation patterns if two leaves arc specimens from two different species (Fig. [I}. Similarly, networks produced 
in silico by optimization routines or developmental simulations that incorporate the effects of biological noise 
exhibit architectures that are to some extent random: each simulation repeat will produce statistically similar, but 
never identical, networks [5Hl2j. To compare naturally occurring networks with the computer simulated models 
we therefore need to be able to test the null hypothesis that the two networks in question have been drawn from 
the same distribution. 

Some of the distribution and structural networks in question are planar, i.e. their edges are (or can be) all 
confined to the same plane and meet only at vertices (no two edges can cross each other). Examples of naturally 
occurring planar networks include the veins of leaves and insect wings, the loopy arterial network of the mammalian 
neocortex and many others. 

Despite the importance of these planar loopy networks, the arsenal of specialized tools and techniques that can 
sufficiently capture the architecture is still limited. Instead, so far the scientific focus has been on quantifying and 
describing the topology of networks with a tree architecture or the connectivity of non-planar complex networks, 
such as the internet. In particular, work developed since the fifties to describe river networks and dendritic 
architectures helped establish some powerful measures to describe the topological properties of tree structures. 
The Horton-Strahler stream ordering system [T31II1] and the Asymmetry [15] are two such measures that played 
a crucial role in understanding the laws that dictate network growth and organization. Invaluable though these 
measures might be for rivers and dendrites, their definition and usage presuposes a tree architecture and loops 
destroy their consistency. 

Although measures developed for general, non-planar complex networks such as the degree distribution and 
the community structure in principle work for planar graphs, frequently they are not fine tuned to capture many 
aspects of the 2-d network organization [TBH22j . Methods to extract the hierarchical organization of complex 
networks have focused primarily on the node connectivity |23j . Similarly, other more specialized metrics such as 
the distribution network entropy [24| are not very informative with regard to quantities of interest in this work, and 
in particular the hierarchical organization of graphs. Some specialized schemes have been developed to quantify 
the loopy architecture of dicotyledonous leaves (see e.g. [25]), and though they can reveal important information 
about leaf physiology and function [26] these methods do not explicitly characterize the nestedness of the topology. 

To achieve a meaningful and elegant quantification of highly interconnected and loopy biological networks we 
need a sufficiently nuanced metric that captures certain important aspects of the topology and architecture of the 
loopy network without relying on the exact value of the bond strength or geometrical location. Such a metric 
would allow phenotypic parameter reduction and assignment of numeric values to the level of loop nestedness and 
other aspects of the architectural organization that are not represented by descriptions that rely on local, scalar 
quantities (such as histograms of the vein density) . More importantly, it will allow a quantitative and topologically 
based comparison between natural loopy networks and the prediction of optimization models. 

In this paper we present a method that allows us to map the architectural organization of a planar graph to that 
of a binary tree. We then use three metrics widely used for binary trees and examine their properties with regard to 
the original graph.. These metrics are the Asymmetry [15], the cumulative size distribution [27]l28] and the Strahler 
bifurcation ratio [T3]. We present results from three classes of networks: computer generated networks (whose 
building rules are predetermined), networks optimized for known functionals and naturally occurring networks such 
as leaf veins and the arterial vasculature of the rat neocortex. We finally discuss the advantages and disadvantages 
of each approach and present future directions and applications. 
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Figure 2: (a) Deletion of an edge in a loopy graph, (i) The deletion of the edge joins two adjacent loops, (ii) The deletion 
of the edge disconnects the graph, (b) Hierarchical decomposition of a planar graph. Boundary loops sequentially join the 
outside space, marked as oo. Left: Nesting tree of the hierarchical decomposition. Right, top to bottom: hierarchically 
decomposed graph. The bottom right panel corresponds to the full graph, the rest represents the network at different 
levels of decomposition (the corresponding cutoff level of the tree representation is marked with a gray dashed arrow). As 
edges of the graph are hierarchically deleted, based on their thickness, the original loops (A-E) are joined to form derived 
loops (N1-N3). (c) Hierarchical decomposition of a planar graph. Phantom boundary loops surround the graph perimeter. 
Loops contiguous to the perimeter of the graph join a ring of phantom boundary loops. The decomposition proceeds as 
in (b), but the phantom loops 61-64 appear among the loops of the original graph in the tree representation, (d) Building 
blocks of a loopy architecture. The two basic building blocks of the loopy architecture can be identified using the tree 
representation of the graph. (ii),(i2): multiplicative nestedness. Nested loops merge hierarchically, (is): This architecture is 
represented by "tall" trees. (iii),(ii2): additive nestedness. Ordered loops join consecutively, (iis): The tree representation 
is that of "short" trees. Graphs (ii) and (12) map equivalently to (is), similarly graphs (iii) and (ii2) map equivalently to 
(lis), (e) Cumulative size distributions of additive and multiplicative models of nestedness. (ii) Nesting tree for additive 
nestedness. The degree of each node is is shown. (12) Degree (size) distribution for additive nestedness. (is) Cumulatize size 
distribution for additive nestedness. (iii) Nesting tree ,(ii2) Degree (size) distribution and (lis) Cumulatize size distribution 
for multiplicative nestedness. 



2 Results 

2.1 Hierarchical decomposition 

We have developed a method that maps a predominately loopy architecture to a dichotomously branching tree. 
This method hierarchically decomposes the loopy architecture by succesively deleting edges and joining contiguous 
loops, and represents this hierarchical decomposition as a tree, termed the nesting tree. 

In what follows, the term link will refer to a graph element that connects two nodes, and the term edge will 
refer to a chain of links, connecting nodes. Each node in an edge is connected to exactly two other nodes, except 
the nodes at the boundaries of the edge, which can be connected to only one other node (when that edge is the 
"leaf" of a tree), or three or more other nodes. The "edge strength" Wj is a quantity that parametrizes the weight 
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of the edge J. If an edge J is composed of a chain of hnks, then Wj can be set to be the edge strength of the 
weakest of the chain hnks, the median value, or any other quantity that is of interest. The term loop is used to 
refer to the graph cycles, and the terminal or ultimate loops are the cycles that do not contain other loops. 

The first step of this hierarchical decomposition framework is a pruning step, where all tree-like components 
rooted on the loopy graph backbone, if present, are removed from the graph. This step eliminates all vertices that 
belong only to one edge, and produces a graph where each edge either separates or connects two loops (Fig. [2][a)). 

In the second step of this hierarchical decomposition, we order the list of graph edges based on their width (if 
the graph edge is composed of a single link) or their edge strength, and identify the edge with the smallest Wj. 
In this step we assume that the edges can be ordered according to their weight in a strictly monotonic fashion, 
namely that Wj ^ Wk for every pair of edges J, K. This is a requirement that can easily be implemented by 
infinitesimally randomly perturbing Wj or Wk when Wj = Wk- 

In the third step, we remove the edge Jg with the smallest edge strength from the graph. When an edge 
separates two contiguous loops, as in Fig. [2Ia)(i), then its removal will result in joining the two loops to form a 
larger one, the area of which is the sum of the areas of the two initial loops. In most cases, this step will also result 
in joining the remaining edges of the contiguous loops. For example, in Fig. [2Ka)(ii), the links AB and BC will be 
joined to form the edge AC. 

We then repeat steps 1,2 and 3 iteratively, to sequentially remove every edge, and as a consequence, gradually 
join every loop, and perform what we have termed a hierarchical decomposition of the graph. We can represent 
this procedure with a dichotomously branching tree, as follows. The "leaves" of the tree are the original loops of 
the full graph, loops A-E in Fig. [Sfb) , and each node downstream of the leaf nodes represents a larger loop that 
is formed by joining two upstream loops through the removal of an edge. The location of the downstream nodes 
on the vertical axis of the branching tree represents the edge strength that was removed to join these two loops. 
Loops are being hierarchically combined until they break to the outer region, termed exterior (and labeled oo). 
The exterior is treated as a separate loop. 

This method will hierarchically decompose the original graph and will register this hierarchical decomposition as 
a binary, nesting tree. The branching patterns of this nesting tree contain information about important topological 
properties of the original graph. The nesting tree allows us to adapt and use measures traditionally used and 
defined on trees, to quantify the architecture and topology of loopy graphs. 

The hierarchical decomposition and the nesting tree contain no explicit information about the geometry of 
each edge and element of the graph, other than the fact that the two joining loops need to be adjacent. Nodes of 
the nesting tree thus correspond to neighborhoods of the original graph - the nesting subtree tj rooted at node j 
represents the architecture of the subgraph enclosed in the loop represented by node j. 

When edges at the graph perimeter are removed and loops at the boundary merge with the exterior, the 
neighborhood information is lost. We can retain that information by appropriately fragmenting the exterior region. 
Instead of having a single exterior loop, where every boundary loop sequentially merges to, we define a multitude 
of exterior loops as follows. We consider an exterior phantom loop that encompasses the original graph in its 
entirety. We then connect the vertices on the perimeter of the original graph with the perimeter of the phantom 
loop as shown on Fig. EJc). Thus defines n boundary phantom loops, labeled &i, 6„, where n is the number of 
loops in the original graph that are adjacent to the perimeter. The added phantom exterior loop and links are 
assigned infinite weights and will never be removed during the hierarchical decomposition. After the addition of 
the phantom loops to the original graph, we proceed to iteratively decompose the graph as before, and represent 
the decomposition with a binary nesting tree, like the one shown in Fig. [2][b). In this way, the neighborhood 
information at the boundaries is preserved and will be reflected in the architecture of the nesting tree. 

The nesting tree facilitates straighforward identification of the two basic building blocks of the organization of a 
planar graph. We will denote these building blocks as multiplicative (Fig. [5][ii)(i2)), and additive (Fig. [2][iii)(ii2)). 
The multiplicative building blocks consist of events where the small loops are joined in an iterative, self-similar 
fashion. It maps to a tall binary tree, such as the one shown in Fig. [IJii ^). The additive building block is 
characterized by sequential joining of minor loops to an encompassing major loop. It maps to a short binary tree, 
as in Fig. [2Kui_2). 

This mapping to a nesting tree is not a bijection. Any information about the geometric organization (shape and 
location) of the loops is lost. Only topological information is retained. For example, networks Fig. [5Jii) and (12) 
both map to Fig.[2ji3), and Fig.[2jiii) and (ii2) both map to Fig.[2jii3). Elements of the architectural organization. 



4 



such as loop area or aspect ratio can be retained by assigning related values to the nodes of the tree j and defining 
quantities that reflect their distribution. For example, the cumulative size distribution is based on measurements 
of the loop areas A{j) assigned to the nodes j of the nesting tree. 

When an edge connects, rather than separates, two loops, its deletion will disconnect the graph (Fig. [2l^a)(ii)). 
There is a number of ways to incorporate such an event to the hierarchical decomposition algorithm. In the 
example cases that we consider in this work, such events are rare, so for simplicity we chose to discard them in our 
implementation. In particular, we replaced the weight value of the disconnecting edges with the maximum edge 
width value of the disconnected loopy components, this way ensuring that the loop will not disconnect from the 
graph before it is hierarchically merged to the encompassing loop (Fig. [2][a)(ii)). 

The nesting tree allows the unique assignment of a number to properties of the hierarchical organization of the 
graph and decouples geometry from topology. In this paper we consider and adapt three measures that have been 
traditionally used on trees: the Asymmetry, the cumulative size distribution and the Strahler bifurcation ratio. We 
apply those measures to the nesting tree and consider what they mean for the organization of the original graph. 



2.2 Asymmetry 

The Asymmetry is a metric that characterizes the topological structure of a binary tree. It was first developed 
mainly in the context of neuronal branching patterns, such as dendritic trees and was defined as the weighted mean 
value of the asymmetry of its partitions. Adjusting the definition and notations of [15| . we define the partition 
asymmetry of a bifurcation vertex j as: 

qir„s,)='-l^ (1) 

with Sj > Vj and Sj + rj > 2. The parameters Vj and Sj are the degrees of the two subtrees at partition j. The 
degree of a (sub)tree is defined here as the total number of the leaf nodes (terminal segments) of that (sub)tree. 
Note that Eq. [T] differs slightly from the definition in [15] . 

The Asymmetry Qritn) of a subtree rooted at node n can now be defined as the weighted average of the 
partition asymmetry q{rj, Sj) of the nodes j £ t„: 

^ d(n)-l 

where j runs over all d{n) — 1 bifurcating vertices of the subtree {d{n) is the degree of the subtree), and Wj is the 
weight of the partition j. Finally, the normalization factor 'w{t„) is defined as: 

d(n)-l 

The averaged Asymmetry Qt{5) of trees of degree 5 is defined as: 

'5tW = ^ E QT{t^) (4) 

"■^ {t,},d(t])=S 

where 715 is the number of nodes with degree 5. In this work, we adjust this definition to be the mean of the 
asymmetry for all the nodes whose degree is within a distance A/2 from 5: 

Qt{5)^ E Q^fe) (5) 

{t,},\d{tj)-5\<A/2 

Calculated on the nesting trees of the hierarchical decomposition, the Asymmetry is an metric that quantifies 
the nestcdness of the original graph. High Asymmetry values correspond to a graph that is primarily composed 
from additive building blocks, and low asymmetry values correspond to a graph that is made from multiplicative 
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building blocks. The actual correspondence between Asymmetry values and level of nestedness depends on the 
choice of weight function Wj . Different choices of weight functions amplify different aspects of the graph architecture, 
and comparisons of Asymmetry plots of different graphs should only be done when the weight function choice is 
consistent. 

In the bulk of this work we will use a weighted averaging window that includes all nodes of the subtree, with 
weight Wj = d(j) — 1. In the Supplemental material we present results acquired by considering no averaging 
Qo{tn) = q{rm Sn) and by averaging over a shallow averaging window. 
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Figure 3: (a) Generated graphs. These graphs were constructed to share identical underlying topology (N = 817 vertices, 
triangular lattice) and edge width distribution, (b) The Asymmetry Qritn) of the every subgraph tn of rooting node n is 
plotted as a function of the base 2 logarithm of the degree d{n), for the nested (red circles), gradient (green squares), and 
peaks model (blue diamonds). For the peaks and random lines model, instances of the graph are plotted with highlighted 
subgraphs of degree 2"^ and 2^ (nested) and ~ 2^ * (peaks). Note the quasi-periodicity of the Asymmetry of the nested 
model (a signature of the self similar structure of the nested model) and the change of monotonicity of the peaks model 
(indicating a qualitative change in the architecture of the graph at that level of organization), (c) The Asymmetry Qritn) 
of the random lines model (red) and random links model (cyan). The x-axis is the logarithm of the degree of the vertex 
or the nesting tree. Red line: averaged Asymmetry of subgraphs of degree d(n), random lines model. Cyan line: averaged 
Asymmetry of subgraphs of degree d{n), random links model. Inset: Density plots: The overlap of the two distributions is 
plotted in white, (d) The averaged Asymmetry Qrid) of the nested (blue), nestedS (orange), nestedlO (light blue), random 
lines (red) and random links model (cyan) as a function of the base 2 logarithm of the degree d. (e) Cumulative size 
distribution P{A > a) of generated models. Random links model (green), nested (blue), gradient (magenta), peaks (green). 
The total area of the graphs has been normalized to 1. Discontinuities or near discontinuities in the slope of cumulative 
size distribution indicate lengthscales where potentially the architectural organization changes qualitatively, (fi). Adjusted 
cumulative size distribution, random links model. (f2)The Adjusted cumulative size distribution P{A > a) * a is plotted 
for the nested (blue), nestedS (orange), nestedlO (light blue) and random lines model (red). The Adjusted cumulative size 
distribution of the self-similar networks (nested, nestedS, and random lines) can be approximated by a straight line of slope 
zero. Notice the periodicity in the nested lines model. The colored area designated the standard error of 20 realizations. 
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2.3 Cumulative size distribution 



The cumulative size distribution j271l28| is the cumulative distribution over the areas associated with the nesting 
tree nodes. It is calculated by assigning an area Value A{j) to each node j of the nesting tree, and then calculating 
the probability P{A > a) that an area drawn at random will exceed a certain value a. In general, we can associate 
the nesting tree nodes with any quantity that reflects a property of the original graph that is of interest, such as 
the total number of terminal loops nested in loop j of the original graph (equal to the degree d{j) of node j of the 
nesting tree, if the terminal loops are of equal size) . 

The cumulative size distribution reflects the overall architecture of the original graph, as the smaller degree 
nodes of an aggressive subdivision, like the one in Fig. [2je) (ii) will be overepresented in the degree and cumulative 
degree distribution. It is easy to show that the cumulative degree distribution of iterative, self similar architectures 
is inversely proportional to the area 

P{A>a)'^l/a. (6) 

Conversely, the cumulative degree distribution of an architecture with additive nestedness (Fig.[2l[e)(i)) is a straight 
line with slope: 

dP{A > a) ^ 1 

da 2' ^ ' 

2.4 Strahler bifurcation ratio 

The Horton-Strahler stream-ordering system has been an invaluable tool in quantifying aspects of river topology 
and architecture since its inception in the fifties by Horton and Strahler [131 [14]. It has since been used with 
considerable success in describing the topology of a wide class of natural and man-made networks. 

According to the Horton-Strahler stream-ordering system, the terminal nodes of the network (the leaves) are 
assigned Strahler order 1. The order of every non-leaf node is determined by the following rule: when two edges 
are connected to two nodes of Strahler order wi,aj2 upstream, the node downstream is assigned an order 

uj = max{uji,Lj2) + Si^-^^uj2- (8) 

The Strahler numbers (or the related Horton numbers) can be used to quantify the tree topology in a number 
of ways. In this work we focus in particular on the Strahler bifurcation ratio, defined as: 

Ruj = Suj/ Suj+i (9) 

where S^j is the number of streams of Strahler order w. A stream is defined as a maximal path of branches 
connecting vertices of Strahler order w, ending in a vertex of higher order. 

The law of stream numbers states that the stream numbers S^) approximate an inverse geometric progression 
with the order w, a statement that implies R^^ = const. However, it is not possible to use this law as evidence of 
self-similarity of a distinctive architecture, as it is followed by the vast majority of binary trees }29) . 

The Horton-Strahler stream-ordering system cannot be directly used to describe loopy networks, as there can 
be no unique assignment of the stream order in a redundant graph. The hierarchical decomposition and the nesting 
tree provide a mapping that allows assignment of Strahler numbers to a loopy graph, as the loops of the original 
graph map to the vertices of the nesting tree and the Strahler number of node j depends on the nestedness of the 
graph segment enclosed by the loop j. 

We now analyze examples from three classes of graphs: models generated by specific, prescribed building rules, 
outputs of optimization routines and natural graphs (in particular the venation of two dicotyledonous leaves and 
the arterial vasculature of the rat neocortex). 

2.5 Hierarchical decomposition of generated networks 

In this section we will consider various classes of architectural models. These computer generated networks were 
produced according to various predetermined rules. We use the hierarchical decomposition and associated metrics 
to quantify various aspects of the architecture, demonstrate what the metrics reveal about the graph organization 
and understand the effects of the finite size, boundaries and of noise. 
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Figure 4: (a) Optimized networks, fluctuations in the load (sink model). Instances of optimized graphs (7 = 0.2,0.5,0.7) 
when the load is concentrated at a single, moving, point, (b) Optimized networks, robustness to damage (bond model). 
Instances of optimized graphs (7 = 0.2, 0.5, 0.7) when robustness is required under the presence of random damage, (c) 
Asymmetry of sink model, (d) Asymmetry of bond model. The average asymmetry Qrid) is plotted as a function of the 
normalized subtree degree d/d^^,^. Red line: 7 — 0.2. Green line: 7 = 0.5. Blue line: 7 = 0.7. Black dashed line: random 
links model. The colored area represents the standard error after averaging over 20 realizations of each model, (e) Adjusted 
cumulative size distribution, sink models. The gray line overlayed on the blue, 7 = 0.7 line is the random links model, (f) 
Adjusted cumulative size distribution, bond models. The adjusted cumulative size distribution P{A > a) * a is plotted for 
7 = 0.2, 0.5, 0.7 (red, green, blue respectively) The adjusted cumulative size distribution is averaged over 20 realizations for 
the bond, sink and random edges model. The colored area represents the standard error after averaging over 20 realizations 
of each model. 



The networks used in this section are presented in FigI3l^a). The underlying geometry, link connectivity and 
point-wise link weight distribution are identical in every model. The architecture is solely defined by the building 
rule according to which the link weight values are assigned on the network. In the gradient model in Fig[3l the link 
weights are distributed according to the link center Euclidean distance from the left-most vertex, creating a smooth 
gradient of link weight. In the peaks model, the thick links are concentrated around seven equidistant peaks. In 
the nested model, the straight lines defined by the underlying link connectivity arc ordered based on a self similar 
subdivision scheme: the lines on the boundaries and center are assigned order k — 1, the lines bisecting order k = 1 
lines arc assigned order k = 2 etc. The link edges are similarly ordered according to weight, and then distributed 
to the ordered straight lines so that higher thickness links occupy lower order lines. This produces a hierarchical 
self-similar pattern, characterized by long range order in the link weights. The nestedd model demonstrates an 
instance of a class of models that is produced by the nested model, choosing five lines at random and randomly 
permuting their order. Similarly, the model nestedlO, is derived from the nested model by swapping 10 lines at 
random. The random lines model is produced by a random permutation of all the lines. Finally, the model random 
links is produced by random assignment of the weights to the links and exhibits no log-range order. 

In Fig. [3]^b) we plot the Asymmetry Qritj) for the architectural models termed nested (blue), gradient (ma- 
genta) and peaks (green). Ignoring the nesting tree vertices that correspond to the phantom boundary loops, the 
nesting tree of the gradient model is a purely additive tree of the type shown in Fig. [Uc)(ii). The nesting tree of 
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the nested model is similar to Fig. [5Je)(iii), however, the iterated building unit is composed of four sequentially 
joining elements, rather than just two joining nodes. 

We have analytically calculated the Asymmetry of Eq. [2] for the infinite gradient and nested models (no bound- 
aries). For the gradient model, it can be trivially found to be: 

QTit,)^l-^y (10) 

For the nested model, the analytical expression in closed form is more complicated and presented in the supplemental 
material. To demonstrate the effects of the boundary in the Asymmetry of the nested and gradient model, we 
overlay the theoretical predictions on the finite size numerical results of Fig.[3Ib). For the gradient model, where a 
large number of low order loops break directly to the boundary, we notice a deviation of the actual measured finite 
size Asymmetry from the theoretical one. This deviation is mostly noticeable for small d. In the nested model 
case, there are no low level loops that join the exterior during the initial stages of the decomposition, so the finite 
size effects produce a deviation from the theoretical graph only at high degrees d. The damped fluctuations in the 
Asymmetry of the nested model are a signature of the model's self similarity. The asymptotic relaxation value of 
these fluctuations is indicative of the iterative building block of this architecture. 

The Asymmetry plot of the peaks model follows closely the one of the gradient model, but, at approximately 
d ~ 2^-^ there is a marked change of monotonicity. This is the characteristic scale where the architecture of the 
model changes qualitatively. Until that scale, the architecture was predominately additive, with smaller loops 
sequentially joining larger ones, and the Asymmetry curve followed qualitatively that of the gradient model. The 
Asymmetry decreases when the six separate, large size segments, represented in the inset graph with different 
colors, join. After those events take place during the hierarchical decomposition smaller loops with stronger edges 
continue to sequentially join creating a pattern in the Asymmetry plot that is again reminiscent of the gradient 
model. The Asymmetry can be used to identify characteristic length scales in graphs where major changes in the 
architecture take place. 

All the three models shown on Fig. EKb) are deterministic, with relatively simple architectures. Models such as 
the random links or the random lines model exhibit a much more complex asymmetry profile, as shown in Fig.[3Kc). 
The asymmetry values in that case are drawn from a distribution the properties of which reflect the architecture 
in question. Calculating mutual information and comparing density maps such as the ones shown in the inset of 
Fig. can provide a statistically meaningful way to examine the null hypothesis if two random graphs belong 
to the same architectural class. An extensive statistical comparison of the different architectural models is beyond 
the scope of this work. Alternatively, we calculate the average Asymmetry plotted in Fig. ^c) with the red 
and cyan solid lines and in Fig. ^d). The exact average asymmetry of each realization of the random models 
depends on the details of the noise. In Fig. ^d) we plot the mean Qrid) over 20 realizations of the nested (blue 
line), nestedS (orange), nestedlO (light blue), random lines (red) and random links (cyan) models. The colored 
area represents the standard error. 

The nestedS and nestedlO models represent intermediate models between the nested and the random lines 
architecture, with progressively increasing disorder (and Asymmetry) as the number of lines that have been swapped 
becomes greater. The nested, nestedS, nestedlO and random lines model are architectures with long range order 
in the link strength, qualitatively significantly different than the random links model, in which the link strength is 
uncorrelated. This difference in reflected in the Asymmetry values of the random lines and random links models 
(Fig-lid)). 

The cumulative size distribution of the generated models is presented in Fig. He) and Fig.dfi 2). In particular, 
in Fig. He) we plot the cumulative size distribution of the peaks (green), gradient (magenta), nested (blue) and 
random links model (cyan). As anticipated, the gradient model follows a straight line of slope 1/2 (a small deviation 
for small a is due to boundary effects). Kinks and discontinuities in the slope, like the ones seen in the peaks model 
curve, are indicative of qualitative changes in the architecture. The random lines and nested model curves are 
significantly different from the gradient model. We can robustly test for scale invariance by defining the adjusted 
cumulative size distribution a ■ P{A > a). According to Eq. [6l for self similar graphs like the nested model we 
expect this quantity to fluctuate around a constant value, and the shape of the fluctuations are indicative of the 
iterative building block of the nested model. 
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Figure 5; (a) Segments of digitized leaf vasculature. The image of the skeletonized leaf has been overlayed with the digitized 
portion of interest, (ai) Bursera teconata, (0,2) Protium heptaphyllum. Images courtesy of Douglas Daly, New York Botanical 
Gardens, (b) Hierarchical decomposition of Bursera and Protium. (ai) Bursera, (a2) Protium. Top to bottom remaining 
loop at three different, progressively higher thickness cutoffs. Notice the persistent minor loops at the proximity of the major 
veins, (c) Segmentation of Protium heptaphyllum and associated tree representation. The protium intercostal area area has 
been separated to six color-coded sectors, as identified by hierarchical decomposition. The associated tree representation 
for that sector is shown for the green and red sector. The non-colored (white) areas of the graph and associated gray links 
on the tree representation correspond to high asymmetry nodes of the tree representation. Note how the high asymmetry 
areas are concentrated near major leaf veins. 



In Fig. EJfi) we plot a ■ P{A > a) for the random links model, and in Fig. ^{2) for the various nested and 
random lines models. As expected, the curves for all realizations of the self similar models fluctuate around a 
straight line. The periodicity of the curve can reveal the size of the architectural unit of the self similar network. 
The deviation from a straight line for large a is due to boundary effects. As the disorder increases, the periodicity 
becomes less pronounced, and disappears at the random lines model. 

2.6 Hierarchical decomposition of optimized networks 

In this section we use the hierarchical decomposition method and the nesting tree to analyze the output of the 
optimization routines presented in [5]. Here, unlike the architectural models presented earlier, the building rules 
according to vifhich the networks were constructed are not a-priori known. However, the functional purpose of the 
networks is known, as they are the (local) minima of global energy functions. The two models under consideration 
are a robustness to damage (broken bond) and fluctuations in the load (sink) model. 

Modeled as electrical (or equivalently water distribution) grids, the networks transport load from the root 
(bottom center vertex in the networks of Fig. Hlja)) to other nodes in the network. In the "bond" model, the 
root has to distribute the load evenly to all the vertices, even if a random single bond is removed (robustness to 
damage). In the fluctuating sink model, instead of a uniform distribution of sinks there is a single sink the position 
of which moves across the network. The cost to build the network is determined by a function K = ^ and is 
set to a constant in each case. The parameter 7 quantifies the "economy of scale" , i.e. how relatively expensive is a 
high conductivity edge compared to a smaller edge. The link thickness of the graphs shown in Fig. \^a) represents 
the bond conductivities, which are determined by optimizing for the total network power dissipation (results are 
shown for 7 = 0.2, 0.5 and 0.7 ). 

The Asymmetry plots demonstrate the strong statistical similarity of the sink 7 = 0.5 and 7 = 0.7 models with 
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the random links model at intermediate and large scales (Fig. |l{c)). For the bond models, the 7 = 0.2 follows 
closely the 7 = 0.5 optimum, and they both exhibit a marked change in monotonicity at larger scales. The overall 
asymmetry increases with 7 in the sink model, whereas there appears to be a significant qualitative change in the 
architecture between the 7 = 0.7 and 7 = 0.2, 0.5 of the bond model. Here it should be noted that the asymmetry 
metric, as defined here, does not depend on the actual numerical value of the bond strengths, just the absolute 
ordering on the lattice. The sink model network for 7 = 0.5, 0.7 appears uniform as the smaller conductivity values 
are similar in value, however, architecturally the network is similar to the random model of Fig. [Sja). 

The adjusted cumulative size distribution shown in Fig|3Jc),(f), overall qualitatively reproduces the findings of 
the Asymmetry. The bond model for 7 = 0.7 exhibits a small size plateau. The sink 7 = 0.7 model follows a similar 
curve as the one of the random links model. Note the change of monotonicity in the bond 7 = 0.2 and 7 = 0.5 
model. This indicates a change of architecture from primarily additive to primarily multiplicative nestedness. 

2.7 Hierarchical decomposition of natural networks 

In this section we apply the methodology we developed so far for two real examples, a leaf from Bursera teconata 
and a leaf from Protium heptaphyllum, show on FigjSl The leaves have been cleared and stained by the group 
of D. Daly in the New York Botanical Gardens, who provided us with high resolution images of the specimens. 
We reconstructed and digitized the vasculature of leaves using custom made software that we have developed to 
translate the pixel values information to a collection of nodes and edges on which we can perform hierarchical 
decomposition. 

In Fig. [5] we show the reconstructed portion of the leaves, overlayed on the digital image from which it was 
acquired. A non-uniform staining or illumination of the specimen can introduce bias to the reconstruction algorithm 
and certain neighborhoods of the reconstructed graph might appear to have spuriously large weights. In particular, 
executing an initial decomposition step on the two networks of in Fig.[5][bi) and (bi), we can easily see that unlike 
the Bursera, the Protium sample appears to have strong loops of smaller size concentrated around major veins. 
A careful inspection of the actual specimen is necessary to determine whether the origin of this bias is due to 
differential staining or this effect is of true biological origin. Although problems like this can be dealt before the 
digitization step in a number of ways (such as a variable threshold), here we will not discuss this, but will accept 
the input data at face value, and discuss cleaning methods that can be applied post the digitization stage. 

A hierarchical decomposition of the intercostal area of Bursera allows us to identify high level nodes of the 
nesting tree that correspond to major loops. We use the nesting tree to identify a natural segmentation of the 
graph to six major areas which we plot in Fig. ^c) along with the corresponding nesting subtrees for two of those 
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Figure 6: (a) Asymmetry of Bursera and Protium intercostal areas. The average asymmetry Qrid) is plotted as a function of 
the normalized subtree degree d. Red solid line: Protium, cleaned. Red dashed line: Protium, full graph. Blue line: Bursera. 
Dark diamonds: random edges model. Dark circles: nested model, (b) Asymmetry of Protium intercostal segments. Qrid) 
is plotted as a function of the normalized subtree degree d. Black dashed line: Protium, cleaned. Red, blue, green, magenta, 
cyan, yellow lines: Protium segments, colorcoded as in Fig. [5] Gray squares: average of segment asymmetry with standard 
error, (c) Adjusted cumulative size distribution, Bursera and Protium. Red solid line: Protium, cleaned. Red dashed line: 
Protium, full graph. Blue line: Bursera. 
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sections. The histogram of the partition asymmetry q defined on the nodes of the nesting tree has a local minimum 
at approximately q ~ 0.97. This value can serve as a natural cutoff for data cleaning, In the nesting trees of 
Fig.ElJc), we color the links of the subtree upstream of the nodes with partition asymmetry higher than 0.97 with 
gray. The corresponding high asymmetry loops are colored white in the original graph. We see that indeed the 
high symmetry loops are consistently concentrated around major veins. 

The asymmetry curve Qt of the intercostal area of Bursera, reaches a plateau. On the contrary, the Protium 
asymmetry does not relax to a constant value. However, if we clean the sample by disregarding the high asymmetry 
nodes with q > 0.97, we see that Protium asymmetry curve similarly reaches a plateau, which is nevertheless 
higher than Bursera, indicating an architectural model based on more additive than multiplicative building blocks 
compared to Bursera. We can calculate the asymmetry for each individual segment of Protium in Fig. [5] and see 
that, as expected, the different segments exhibit the same architecture and the asymmetry curves relax to a value 
of approximately Qrid — !> 1) ~ 0.6, significantly different than the value of 0.45 of the Bursera. 

The cumulative size distributions of Fig. [5] qualitatively follow the observations from the asymmetry plots. The 
cleaned Protium curve, as well as the Bursera curve, both reach a plateau, however the cumulative asymmetry 
cannot effectively distinguish between the two speciments. 
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Figure 7: Strahler bifurcation ratio for the various generated, optimized and natural graphs presented in this paper. Red 
error bar: standard error of the linear regression fit (represents goodness of linear fit). Black error bar: standard deviation 
of the logarithm of the bifurcation ratio (average over 20 realizations). Insets: Number of Strahler streams Su of order uj as 
a function of uj for the random lines, nested and gradient model and the Bursera leaf. Note that in each case, the 5*^ follows 
closely an inverse geometric progression with to (shown with the red dashed line). 
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2.8 Strahler bifurcation ratio 



The Strahler bifurcation ratio (jH]) defined on the nesting tree can provide a metric to quantify the overaU nestedness 
of graphs. Since the Strahler law of stream numbers is an inevitable reality for most trees, it is possible to fit the 
plots log(S'tj) versus logoj with a straight line the slope of which will determine the logarithm of the Strahler 
bifurcation ratio Rs for the whole graph. Examples of this fit are shown in the inset of Fig. [71 The best fit is found 
in the least squares sense, and it is forced to pass through (1, S'l) (S'l is equal to the total number of ultimate loops, 
or leaf nodes in the nesting tree). The data point for max(cj) is discarded, as it is very sensitive to noise. 

In Fig. [7] we plot the Strahler bifurcation ratios for all the graphs presented in this paper. For the architecture 
or the optimization models that are not deterministic each realization of the graph will produce a different Rs- 
In those cases, we plot (Rs), the average bifurcation ratio over 20 realizations, with the black error bar being the 
standard deviation. The red error bar represent the (average) goodness of the linear fit. Notice the extent of the 
red error bar for the gradient and peaks models. 

The Strahler bifurcation ratio can clearly distinguish between the strongly multiplicatively nested Bursera and 
additively nested gradient model, but, with our current implementation it could not sufficiently distinguish between 
many of the models presented in this work. A major drawback of Rs is that it is a single number which is inherently 
unsufficicnt to capture the complexity of networks whose architectural properties do not necessarily remain the 
same over all lengthscales. 

2.9 The rat brain 

The analysis and methodology presented in this work can be useful not only for leaves, but any other, biological 
or man made, planar graphs. A notable example is the arterial vasculature of the rodent neocortex which forms a 
planar network with multiple loops [3] . We extracted the diameters of the arterial blood vessels from a composite 
rat brain image provided to us by the Kleinfeld group in UCSD and augmented the connectivity information in [3] 
to obtain a weighted map of the arterial vasculature of the rat brain. Although the resolution of the image in our 
disposal does not allow us to determine the vein widths with absolute confidence, we were able to identify major 
vascular sectors and determine that, according to the data at hand, the architecture of the network in question is 
primarily additive than multiplicative. 

3 Discussion 

We have presented a framework that allows us to quantify the hierarchical organization of predominately loopy 
architectures. Our hierarchical decomposition consists of three iteratively repeated steps: 

1. pruning of the tree-like components 

2. ordering of the edges 

3. removal of the thinnest edge 

This framework relies on the mapping of loopy planar graphs and their hierarchical decomposition to binary nesting 
trees. The nesting tree is subsequently used to quantify the architectural organization of the original graph. A 
number of quantities that reflect various aspects of the graph organization can be defined on the nesting tree, 
each with each own advantages and disadvantages. In this work we presented results for three such quantities, the 
Asymmetry Qt-, the cumulative size distribution and the Strahler bifurcation ratio. The Asymmetry is a bottom- 
up approach that assigns a number to every composite loop at each scale. This number is a weighted average of 
the nestedness of the architecture of the portion of the graph enclosed in that loop. Depending on the averaging 
window, two different architectures of a high degree loop can map to the same Qt value. On the contrary, the 
cumulative size distribution performs better in differentiating architectures at the high levels of organization. The 
larger number of low level loops frequently results in washing out interesting features of the structures at smaller 
scales. 

These observations are demonstrated in the sink and bond model Asymmetries and cumulative distributions of 
Fig. m For example, the Asymmetry of the sink 7 = 0.7 and bond models (Fig. lUJc), (d)) has a local maximum. 
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a feature that is absent from the adjusted cumulative size distribution (Fig. Ul^e), (f)). Similarly, the Asymmetry 
of all bond models is indistinguishable for large scales, whereas the cumulative size distribution can statistically 
distinguish these models. 

Depending on the weight function Wj , the Asymmetry can be used to define a single number that encompasses 
information about the whole architectural organization (e.g. by calculating Qrid — >■ dmax)) rather than to examine 
the architecture at all levels of organization. Such a number would be meaningful only for graphs with some 
degree of self-similarity. The Strahler bifurcation ratio Rg is is used to describe the overall architecture, but it does 
not perform well for complex architectures. We have examined the Strahler bifurcation ratio as a function of the 
Strahler order uj and degree d as an attempt to extract information about the scale dependent organization of the 
graph. We have found that the result is very sensitive to noise, especially at high w. 

The metrics presented in this paper focus on the metric topology of the structure but they do not explicitly 
capture any information about the geometry of the network. The cumulative area distribution depends on the 
area of the terminal loops (the areoles of a leaf vein network). The cumulative area distribution follows closely the 
cumulative degree distribution provided that the terminal loops are not substantially polydisperse. It is evident 
we can supplement the descriptions presented here with more detailed geometrical analysis, in which some aspects 
of the geometry of the closed loops is kept, such as e.g. an approximating SVD ellipsoid, which can be used to 
define a major axis and an excentricity. We can then incorporate such geometrical information into the analysis 
of nesting, i.e., relationships defining what is the average orientation of subloops in relation to the parent loops. 
Such detailed geometrical analysis, however, will evidently be subordinate to the coarser topological analysis we 
have presented here. 

A big part of our extensive understanding of fluvial networks is due to the development of methodology to 
characterize and quantify tree architectures. Accordingly, progress in understanding loopy networks, which are 
ubiquitous in both natural and man made structures, is contingent on our ability to measure their hierarchical 
architecture. The approach presented in this work provides a robust mathematical description of the network 
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Figure 8: Digitized arterial vasculature of rat neocortex and corresponding nesting tree representation. (a)The arterial 
network forms a planar graph. Different segments of the network, as identified by hierarchical decomposition are represented 
by different colors, (b) Nesting tree of the digitized network, the highlighted segments of the network are color-coded. 
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architecture, applicable to leaf venation and other loopy distribution (and structural) structures. It can be used to 
characterize the in silico networks obtained from computer simulations as well as to perform quantitative statistical 
comparisons between theory and experiment. As such, it can provide an invaluable tool in deciphering the functional 
signficance of the loopy networks and possibly their developmental origin. 
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